/*
- 
- Yuri Toledo Neves 8124361
- Trabalho de Recuperação SSC0143 - Programação Concorrente 
- Versão Sequencial
- 
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>


//Funcoes de ajuste
double forcaArrasto (double H, double v);
double ajusteMassa (double t);
double ajusteForca (double t);
double ajusteGravidade (double H);

//EDOs 
double funcaoAceleracao(double t, double v, double gama, double H);
double funcaoRotacao(double t, double v, double gama, double H);
double funcaoHorizontal(double v, double gama);
double funcaoAltitude(double v, double gama);

const double h = 0.0001;	// h do método de Runge-Kutta
const double g0 = 9.80665;	// Gravidade
const double R = 6335437.00000;	// Raio da terra
const double ro0 = 1.2922;	// Densidade do ar no nível do solo
const double hyp = -0.000101741; 	// Resultado equacao hypsométrica
const double PI = 3.141592;	//PI 

double	m0, mp, tq, fempuxo, alpha, coef_arrasto, sf, hf;	//Váriaveis de entrada
double 	r = 0.01; // Raio do foguete

int main(int argc, char const *argv[])
{
	
	double erro;	// Erro
	int i;
	double altura = 0;
	clock_t start, end; 

	double gama_graus = 65, gama = 0;	// Angulo de lançamento	
	double H = 0; 	// Altura 
	double v = 0.1; // Velocidade
	double t = 0;	// Tempo
	double S = 0;	// Deslocamento vertical 

	double k1, k2, k3, k4; // Coeficientes do método para funcao Aceleracao
	double l1, l2, l3, l4; // Coeficientes do método para funcao Rotacao
	double m1, m2, m3, m4; // Coeficientes do método para funcao Horizontal
	double q1, q2, q3, q4; // Coeficientes do método para funcao Altitude 

	if(strcmp(argv[1], "G") == 0)
	{	
		m0 = atof(argv[2]);
		mp = atof(argv[3]); 
		tq = atof(argv[4]);
		fempuxo = atof(argv[5]);
		alpha = atof(argv[6]);
		coef_arrasto = atof(argv[7]);
		r = atof(argv[8]);
		sf = atof(argv[9]);
	}

	if(strcmp(argv[1], "R") == 0)
	{	
		m0 = atof(argv[2]);
		mp = atof(argv[3]); 
		tq = atof(argv[4]);
		fempuxo = atof(argv[5]);
		alpha = atof(argv[6]);
		coef_arrasto = atof(argv[7]);
		sf = atof(argv[8]);
	}

	start = clock();	// Inicio da computação dos resultados

	do{
		
		if(strcmp(argv[1], "R") == 0)
		{	
			if(gama_graus > 70)
			{
				gama_graus = 65;
				r += 0.01;
			}
		}

		altura = 0;
		v = 0.1; // Velocidade
		t = 0;	// Tempo
		S = 0;	// Deslocamento vertical 
		H = 0;
		gama = (gama_graus*PI)/180;	// Converte o angulo de graus para radianos 
		
		// Runge-Kutta de Ordem 4
		for(i = 0; altura >= 0 ; i++)
		{
			
			k1 = funcaoAceleracao(t, v, gama, H);
			l1 = funcaoRotacao(t, v, gama, H);
			m1 = funcaoHorizontal(v, gama);
			q1 = funcaoAltitude(v, gama);

			k2 = funcaoAceleracao(t+(h/2), v+(h/2)*k1, gama+(h/2)*l1, H+(h/2)*q1);
			l2 = funcaoRotacao(t+(h/2), v+(h/2)*k1, gama+(h/2)*l1, H+(h/2)*q1);
			m2 = funcaoHorizontal(v+(h/2)*k1, gama+(h/2)*l1);
			q2 = funcaoAltitude(v+(h/2)*k1, gama+(h/2)*l1);

			k3 = funcaoAceleracao(t+(h/2), v+(h/2)*k2, gama+(h/2)*l2, H+(h/2)*q2);
			l3 = funcaoRotacao(t+(h/2), v+(h/2)*k2, gama+(h/2)*l2, H+(h/2)*q2);
			m3 = funcaoHorizontal(v+(h/2)*k2, gama+(h/2)*l2);
			q3 = funcaoAltitude(v+(h/2)*k2, gama+(h/2)*l2);

			k4 = funcaoAceleracao(t+h, v+(h*k3), gama+(h*l3), H+(h*q3));
			l4 = funcaoRotacao(t+h, v+(h*k3), gama+(h*l3), H+(h*q3));
			m4 = funcaoHorizontal(v+(h*k3), gama+(h*l3));
			q4 = funcaoAltitude(v+(h*k3), gama+(h*l3));

			t = t + h;
			v = v + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
			gama = gama + (h/6)*(l1 + 2*l2 + 2*l3 + l4);
			S = S + (h/6)*(m1 + 2*m2 + 2*m3 + m4);
			H = H + (h/6)*(q1 + 2*q2 + 2*q3 + q4);

			altura = H;
		}

		gama_graus += 0.1;

	}while(fabs((S - sf)/sf) > 0.01);
	
	end = clock();	// Fim da busca pelos resultados

	erro = fabs((S - sf)/sf);
  printf("%s %.5lf %.5lf\n", argv[1], (gama_graus-0.1)*PI/180, r);

	return 0;
}

double forcaArrasto (double H, double v)
{

	double ro = ro0*exp(hyp*H);		// Densidade
	double D  = 0.5*(coef_arrasto*PI*r*r*ro*v*v);
	
	return D;

}

double ajusteMassa (double t)
{
	if(t < tq)
		return(m0 - mp*(t/tq));
	else
		return(m0 - mp);
}

double ajusteForca (double t)
{
	if(t <= tq)
		return(fempuxo);
	else
		return(0);
}

double ajusteGravidade (double H)
{
	return (R*R*g0)/((R+H)*(R+H));
}

//EDO da Aceleração
double funcaoAceleracao(double t, double v, double gama, double H)
{
	double arrasto = forcaArrasto(H, v);
	double massa = ajusteMassa(t);
	double aceleracao = ((ajusteForca(t)*cos(alpha))/massa) - (arrasto/massa) - (ajusteGravidade(H)*sin(gama));

	return aceleracao;
}

//EDO da Rotação
double funcaoRotacao(double t, double v, double gama, double H)
{
	double rotacao = (((ajusteForca(t)*sin(alpha))/(ajusteMassa(t)*v)) - ajusteGravidade(H))*(cos(gama)/v);
	return rotacao;
}

//EDO deslocamento Horizontal
double funcaoHorizontal(double v, double gama)
{
	double horizontal =  v*cos(gama);
	return horizontal;
}

//EDO Altitude
double funcaoAltitude(double v, double gama)
{
	double altitude = v*sin(gama);
	return altitude;
}


	